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INTRODUCTION 


Boundary  element  methods  (BEM)  have  become  an  accepted  alternative 
to  domain-based  methods  (finite  element  and  finite  difference)  for  many 
classes  of  boundary  value  problems.  The  direct  and  indirect  BEM 
formulations  are  two  reasonably  established  formulations  that  are  used 
to  solve  a  variety  of  boundary  value  problems  in  engineering.  These 
methods  require  only  the  boundary  to  be  subdivided  and  are  inherently 
suitable  when  modeling  infinite  domains.  Since  the  methods  are 
relatively  new,  they  have  not  been  as  extensively  developed  as  domain 
methods  for  nonlinear  applications.  They  also  lack  the  generality  (in 
terms  of  extensive  continuum  and  structural  element  libraries)  that 
commercial  finite  element  method  (FEM)  computer  programs  possess. 

Considering  the  individual  strengths  of  domain-  and  boundary-based 
methods,  some  classes  of  problems  (e.g.,  nonlinear  soil-structure  inter¬ 
action)  may  be  most  effectively  solved  by  combining  the  methods. 

Objective 

The  ultimate  objective  of  this  research  is  to  determine  whether  or 
not  the  boundary  element  and  finite  element  solution  methods  can  be 
combined  with  advantage  towards  economical  solutions  of  nonlinear  struc¬ 
tural/geotechnical  problems.  However,  our  immediate  emphasis  is  on 
improving  the  accuracy  of  the  indirect  boundary  element  method.  The 
accuracy  of  a  quadratic  isoparametric  boundary  element  with  a  new 
adaptive  integration  technique  is  evaluated. 

Background 

An  early  application  of  the  direct  method  for  elasticity  is  given 
by  Cruse  and  Rizzo  (Ref  1).  In  this  formulation  Betti's  theorem  is  used 
to  transform  a  volume  integral  to  a  surface  integral.  The  unknown 


boundary  values  for  traction  and  displacements  are  solved  for  directly 
in  the  system  of  equations.  Internal  responses  (stresses,  strains,  or 
displacements)  are  obtained  as  integrals  of  the  tractions  and  displace¬ 
ments  on  the  boundary  of  the  domain. 

An  early  application  of  the  indirect  method  for  elasticity  is  given 
by  Massonnet  (Ref  2).  In  the  indirect  formulation  the  actual  boundary 
value  problem  is  replaced  by  an  auxiliary  problem  in  the  infinite  plane 
(two-dimensional),  and  artificial  boundary  tractions  that  satisfy  the 
boundary  conditions  and  the  governing  differential  equations  constituting 
the  actual  problem  are  solved.  These  artificial  tractions  are  then 
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integrated  to  obtain  internal  and  boundary  responses  of  the  actual 
problem. 

Early  efforts  used  analytically  integrated  elements  with  constant 
or  linear  interpolation  of  boundary  variables.  Lachat  and  Watson  (Ref  3) 
incorporated  isoparametric  element  representations  typical  of  the  finite 
element  method  (Ref  4).  Geometry  as  well  as  boundary  values  could  be 
interpolated  at  a  consistent  order  (linear,  quadratic,  or  cubic)  over  a 
given  element. 

Isoparametric  representations  prevent  closed-form  integration  over 
individual  elements  and  numerical  integration  is  required.  These  element 
integrations  represent  the  main  computational  effort  of  the  BEM.  Thus, 
there  has  been  much  effort  aimed  at  efficient  numerical  integration 
(Ref  5,  6,  and  7).  Lachat  and  Watson  (Ref  5)  present  a  variable-order 
numerical  integration  strategy  with  provisions  for  element  subdivision 
when  the  required  order  of  integration,  for  a  prescribed  accuracy, 
exceeds  that  available  in  the  program. 

A  Naval  Civil  Engineering  Laboratory  (NCEL)  study  of  the  BEM  (Ref  8) 
compared  the  direct  and  indirect  methods  in  two  dimensional  elastostatics 
with  constant  distribution  elements.  The  direct  method  performed  poorly 
(as  compared  to  the  indirect  method)  in  a  region  near  the  boundary  but 
gave  good  values  on  the  boundary.  This  poor  performance  in  the  near 
boundary  region  was  due  to  two  factors:  (1)  the  direct  method  computer 
code  calculated  the  element  integrations  with  four-point  Gauss  quadrature 
as  compared  to  the  indirect  method  computer  code,  which  contained  closedform 
integration  formulas;  and  (2)  internal  response  calculations  in  the 
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direct  method  must  integrate  a  function  with  a  higher  order  singularity. 
Based  on  the  performance  comparison  of  that  study  and  the  fact  that  the 
indirect  method  is  less  formal  theoretically  and  simpler  to  physically 
comprehend,  the  indirect  method  was  pursued  further. 

A  following  study  (Ref  9)  investigated  a  preliminary  coupling  of 
the  finite  element  and  boundary  element  methods  in  two-dimensional 
elastostatics .  The  methods  were  coupled  within  the  computational  frame¬ 
work  of  the  BEM,  satisfying  compatibility  and  equilibrium  explicitly  at 
the  domain  interfaces.  Qualitatively  the  coupling  was  a  success,  but 
the  simple  constant  elements  used  for  both  methods  prevented  substantial 
quantitative  evaluation.  Coupling  the  two  methods  conceptually  as  two 
BEM  domains  provided  a  good  initial  study  but  was  not  as  applicable  to 
the  long-term  goals  for  application  in  nonlinear  soil-structure  interaction 
because  the  approach  yielded  a  system  of  unsymnetric  algebraic  equations. 
Two  areas,  therefore,  required  further  investigation:  (1)  accuracy  of 
the  BEM  formulation,  and  (2)  a  FEM-BEM  coupling  within  a  FEM  computational 
framework. 

Scope 

The  present  study  addresses  the  accuracy  of  the  indirect  BEM.  The 
main  focus  of  the  study  is  a  new  algorithm  based  on  a  recursive  subdivision 
of  boundary  elements  that  improves  the  accuracy  of  isoparametric  element 
integrations  in  the  near-boundary  region.  Numerical  integration  error 
is  studied  in  the  context  of  the  error  field  near  a  single  boundary 
element  in  an  infinite  plane.  Based  on  this  numerical  study,  a  criterion 
for  element  subdivision  is  determined  as  a  function  of  desired  accuracy. 

The  recursive  algorithm  is  presented  along  with  some  initial  numerical 
examples  to  illustrate  its  effectiveness.  Also  included  are  the  refinement 
of  a  variable-order  integration  formula,  which  is  commonly  used  in 
boundary  element  work,  and  an  investigation  of  numerical  error  that  is 
characteristic  of  the  indirect  BEM  near  boundary  discontinuities. 

In  addition  to  numerical  integration  error,  the  indirect  boundary 
element  method  possesses  error  that  is  traceable  to  its  inherent  use  of 
artificial  boundary  tractions.  The  actual  boundary  responses  (tractions 
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or  displacenents)  may  have  small  gradients,  while  the  corresponding 
artificial  boundary  tractions  may  have  very  high  gradients.  Thus,  a 
boundary  subdivision  that  appears  sufficiently  adequate  for  the  boundary 
response  can  fail  to  accurately  model  the  artificial  tractions  and 
thereby  result  in  poor  accuracy  for  internal  response  calculations. 

These  errors  occur  near  geometry  and  traction  boundary  discontinuities. 
One  numerical  example  presented  includes  a  study  of  this  error  as  well 
as  the  error  from  numerical  integration. 

The  computational  aspects  of  this  study  were  carried  out  on  a 
microcomputer.  The  main  incentives  for  this  approach  were  limitation  of 
cost  and  the  availability  at  NCEL  of  an  effective  software  development 
environment  (the  University  of  California  at  San  Diego  (UCSD)  p-System) 
and  a  modern  language  (UCSD  Pascal).  Our  use  of  a  modern  language  has 
led  us  to  pursue  approaches  such  as  the  use  of  recursive  constructs  that 
would  not  have  been  considered  in  a  FORTRAN  environment.  Strongly 
data-typed  languages  such  as  Pascal  can  generally  provide  a  good  research 
environment  because  of  their  enhanced  data  structures,  internal  documen¬ 
tation,  and  modularity. 

Our  recursive  algorithm  is  implemented  for  the  indirect  BEM  and 
with  a  quadratic  isoparametric  boundary  element.  The  concepts  underlying 
the  indirect  formulation  are  reviewed  in  the  following  section. 


INTEGRAL  EQUATIONS 

The  applicable  integral  equations  are  given  here  for  reference  and 
perspective.  Detailed  accounts  of  the  integral  equation  development  are 
given  by  Banerjee  and  Butterfield  (Ref  10)  and  Crouch  and  Starfield 
(Ref  11).  The  equations  given  in  the  following  section  closely  follow 
those  given  by  Banerjee  and  Butterfield  (Ref  10).  Body  forces  and  rigid 
body  translations  are  omitted  for  brevity. 

Both  the  direct  BEM  and  indirect  BEM  are  formulated  in  terms  of  the 
fundamental  singular  solution,  the  Kelvin  solution  for  linear,  isotropic, 
plane  strain  conditions,  which  expresses  the  displacement  field  u^(x) 
due  to  a  unit  force  e^O-  The  indices  i,  j,  and  k  assume  values  of  1 
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or  2,  and  repeated  indices  imply  summation  in  what  follows.  The  Cartesian 
coordinates  x  and  £  represent  the  field  and  source  points,  respectively. 
The  Kelvin  solution  is  given  by 
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If 


u.(x)  =  G.k(x,o  ek(4) 
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A^k  =  arbitrary  constant  tensor  based  on  zero  displacement 
reference  distance 

y .  =  x.  -  £. 


6..  =  Kronecker  delta  function 
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By  incorporating  the  linear  strain-displacement  relationship  and 


Hooke's  law,  the  stress  field  a„(x)  is  given  as 


aij(x)  =  Tijk(x’C)  ek(^ 


(2a) 


where 
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Equilibrium  conditions  applied  at  a  boundary  point,  indicated  by  a 
unit  outward  normal  n^vx),  and  Equation  2a  combine  to  give  the  surface 
tractions  t^ (x)  as 
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|  t.00  =  F.k(x,C)  ek(|) 

l  where 
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(3b) 


Figures  1  and  2  illustrate  the  singular  behavior  of  the  fundamental 
solutions  G  and  T  ,  respectively,  for  a  point  load  applied  at  the 
origin  of  the  Cartesian  system.  The  singularity  of  G  is  order  ln(r), 
while  the  singularity  of  T,  obtained  from  derivatives  of  G,  is  order 
1/r. 

Consider  the  elastostatic  boundary  value  problem  shown  in  Figure  3 
for  a  linear,  isotropic,  homogeneous  domain  Q  subjected  to  the  traction 
and  displacement  boundary  conditions: 


t.(x) 

=  t.(x) 

on  r 

(4a) 

U.(x) 

=  U.(x) 

on  F 

u 

(4b) 

where  t^(x)  and  u^(x)  are  prescribed  distributions  of  boundary  tractions 
and  displacements,  respectively.  The  boundary  conditions  are  not  mutu¬ 
ally  exclusive;  any  two  of  the  four  boundary  values  are  prescribed  in  a 
piecewise  continuous  manner  for  a  well-posed  problem.  Boundary  values 
can  be  defined  in  any  orthogonal  coordinate  system  and  are  not  necessarily 
limited  to  the  global  x^x^  system. 

When  using  the  indirect  BEM,  the  domain  Q  is  embedded  in  an  infinite 
plane  as  shown  in  Figure  3.*  Artificial  tractions,  Pk(£),  acting  on  the 
boundary  F  are  sought  that  satisfy  the  prescribed  boundary  conditions. 

The  artificial  tractions  can  be  expressed  as  a  "continuous"  distribution 


*The  plane  strain  formulation  can  be  converted  to  the  plane  stress 
formulation  by  specifying  an  effective  Poisson  ratio  V  =  v/(l  +  v) . 
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of  unknown  point  loads  analogous  to  e^CO  and  related  to  field  variables 
by  Equations  1  through  3.  The  displacements  and  stresses  for  any  point 
x  due  to  the  artificial  tractions  are  then  given  by  integrals  over  the 
boundary 


Ui(x)  =  l  Gik(x’£}  Pk(€) 

Vx)  =  -T  Tijk(x^)  pk(C)  dC 


a  e  D 


The  governing  differential  equations  are  satisfied  over  the  entire  plane 
including  the  domain  fl,  since  the  responses  are  expressed  as  a  super¬ 
position  of  the  fundamental  solution.  From  Equations  3a  and  3b  the 
tractions  acting  on  a  tangent  line  defined  by  a  normal  vector  n.  are 
given  by 

t.oo  =  f  f..(x,4)  p. a)  di  (7) 


If  the  field  point  x  is  now  allowed  to  approach  the  boundary  T  and  the 
boundary  conditions,  Equations  4a  and  4b,  are  enforced,  we  obtain  boundary 
integral  equations  relating  the  known  boundary  values  in  terms  of  the 
unknown  artificial  boundary  tractions  as  follows, 


ti(x)  =  j.  Fik(x’£}  Pk(^  °n  rt 


Ui(x)  =  Gik(x’^  Pk(^  ^  °n  ru 


Equation  8  must  be  interpreted  as  a  Cauchy  principal  value  integral  and 
is  thus  written  as 


assuming  a  tangent  line  through  x. 


NUMERICAL  FORMULATION 


Equations  9  and  10  represent  an  exact  boundary  integral  equation 
formulation  of  the  problem  of  determining  artificial  tractions.  Practical 
engineering  methods  require  that  the  equations  be  solved  numerically. 

Two  numerical  approximations  are  made: 

(1)  The  boundary  integrations  are  performed  in  a  piecewise  manner 
over  discrete  subdivisions  of  the  boundary  called  boundary 
elements.  Within  each  element  tractions  (real  and  artificial), 
displacements,  and  geometry  are  usually  interpolated  by  poly¬ 
nomial  functions. 

(2)  Equations  9  and  10  are  satisfied  at  discrete  points  within 
each  element  by  the  collocation  method.  This  can  be  considered 
as  a  weighted  residual  method  with  Dirac  delta  weighting 
functions  having  origins  placed  at  the  collocation  points, 
i.e.,  points  where  the  equations  are  exactly  satisfied. 

The  integral  equations  can  thus  be  approximated  by  a  set  of  linear, 
simultaneous  algebraic  equations. 

First  consider  the  piecewise  integration  over  the  boundary.  If  the 
boundary  is  divided  into  Q  boundary  elements,  the  displacement  at  a 
response  point  x  is  given  as 

Q 

^(x)  =  I  J  G.k(x,C)  Pk(£)  d£  (11) 

^=1  ATq 

In  an  isoparametric  element  formulation  any  field  variable  y  (e.g., 
representing  boundary  geometry,  traction,  or  displacement)  defined  on  an 
element  q  can  be  interpolated  in  terms  of  nodal  values  and  shape  functions 
(Ref  5)  as  follows 


Yi  = 
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or=l 


n  (n)  y. 

a  '  'i a 
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(12) 


where  or  is  the  node  index,  M  is  the  number  of  nodes  per  element,  rj  is 
the  normalized  (-1  <  r)  <  1)  local  curvilinear  coordinate,  y^a  are  nodal 
values  of  the  field  variable,  and  (rj)  are  the  appropriate  polynomial 
shape  functions  for  the  element.  A  quadratic  element  (M=3)  is  shown  in 
Figure  A.  The  shape  functions  for  this  element  are  defined  as 


N1 

=  y  (n2  -  n) 

(13a) 

N2 

=  a  -  n2) 

(13b) 

N3 

=  y  (n2  +  n) 

(13c) 

By  interpolating  the  distribution  of  artificial  tractions  on  element  q, 
P^,  with  Equation  12,  the  displacement  at  point  x  can  be  written  as* 

Q  M  1 

u  .00  =  21  P?  /  G?  N  Jq  dn  (14) 

1  q=l  or=l  K"  -1  1K  " 


where  J^  is  the  Jacobian  relating  the  global  Cartesian  system  £  to  the 
local  curvilinear  system  r|-  The  Jacobian  is  given  by 


It  is  a  function  of  rj ,  and  its  value  is  an  indicator  of  the  geometric 
distortion  of  the  qth  boundary  element.  Introduction  of  the  shape 
functions  allows  the  unknown  nodal  values  of  artificial  traction  Pjj  tc 
be  removed  from  the  integrand  as  constant  coefficients  of  the  integral. 

An  expression  for  the  stress  and  traction  at  any  point  x  can  be 
written  in  a  similar  manner  as  follows 


*The  superscript  q  is  not  to  be  interpreted  as  an  exponent . 


q=l  a=l 


?2  /  Fq  N  Jq  dn 

ka  lk  or  1 


Equations  14  and  17  can  be  applied  at  collocation  points  on  the  boundary 
to  obtain  an  approximate  solution  to  Equations  9  and  10.  For  isoparametric 
elements  the  number  of  element  collocation  points  is  equal  to  the  number 
of  nodal  points.  The  collocation  points  can  be  positioned  at  the  geometric 
node  points  giving  continuous  elements,  or  the  points  can  be  positioned 
within  the  element  giving  discontinuous  or  nonconforming  elements  (Ref  12). 
The  known  boundary  values  are  interpolated  at  the  collocation  points, 
thus  giving  a  system  of  algebraic  equations  in  terms  of  the  unknown 
tractions  Pq  .  After  having  solved  this  system  for  Pqa,  the  displacement 
and  stress  response  in  the  domain  Q  or  on  the  boundary  T  can  be  obtained 
by  applying  Equations  14  and  16,  respectively. 

The  element  integrations  of  Equations  14,  16,  and  17  constitute  the 
main  computational  effort  of  the  indirect  boundary  element  method  for 
linear  elastostatic  problems  and  are  the  focus  of  this  study. 


ELEMENT  INTEGRATION 

Except  when  using  simple  elements,  numerical  integration  is  generally 
a  necessary  part  of  the  procedure  for  calculating  the  system  of  equations 
and  internal  responses  in  the  BEM.  In  developing  the  system  of  equations, 
when  the  collocation  point  is  on  the  boundary  element  to  be  integrated 
the  integration  must  receive  special  treatment.  These  integrations  over 
the  singularity  are  treated  in  detail  by  Watson  (Ref  13)  and  by  Banerjee 
and  Butterfield  (Ref  10).  A  combination  of  analytical  and  numerical 
integration  is  usually  required. 


'•rYs.'+y/s/: 
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The  element  integrations  for  internal  response  calculations  can 
also  require  special  treatment.  We  will  concentrate  on  the  evaluation 
of  stress  because  it  is  associated  with  a  stronger  singularity  than 
displacement  (compare  Figures  1  and  2).  One  approach  to  the  integration 
is  the  single  application  of  Gaussian-Legendre  numerical  integration 
(Ref  14).  For  a  given  function  k(r|),  its  integral  is  numerically  evalu¬ 
ated  as 

1  M 

/  k(n)  dn  =  I  A.  k(n.)  (18) 

-1  i=l  1 

where  and  r)^  are  the  integration  weights  and  positions  (Gauss  points), 
respectively,  for  integration  order  M.  In  the  indirect  BEM  this  represents 
the  replacement  of  a  continuous  traction  with  a  distribution  of  discrete 
forces  at  positions  rj^.  The  stress  field  is  thus  given  by  the  superposition 
of  singular  fields  of  the  type  shown  in  Figure  2  with  origin  placed  at 
each  integration  point  and  with  corresponding  weight  factor  A^ 

The  accuracy  of  element  integrations  can  be  examined  experimentally 
in  the  context  of  a  single  element  in  the  infinite  plane,  as  shown  in 
Figure  5.  Figure  5  also  depicts  the  response  region  and  view  direction 
for  graphical  results  presented  in  Figures  6,  7,  and  8.  Figures  6a,  7a, 
and  8a  show  the  stress  (a^)  field  in  the  response  region  using  various 
orders  of  numerical  integration.  Four-point  integration  is  typically 
used  in  codes  where  the  order  of  integration  is  not  designed  to  vary. 

The  singular  behavior  at  the  Gauss  points  is  clearly  evident.  As  the 
element  (which  lies  along  the  line  x=0)  is  approached  the  stress  and 
thus  the  error  become  unbounded  near  a  Gauss  point.  The  error  in  numerical 
integration  can  be  evaluated  by  comparing  the  results  with  those  obtained 
from  an  analytically  integrated  constant  element.  The  cr^  error  field 
for  various  orders  of  integration  is  shown  in  Figures  6b,  7b,  and  8b.  A 
positive  value  indicates  an  overestimate  of  the  compressive  stress 
magnitude.  The  error  plots  are  arbitrarily  clipped  at  ±2%.  The  variation 
in  weight  factors,  A^,  for  the  integration  points  are  reflected  in  the 
variation  in  the  peaks  on  all  the  figures.  At  half  an  element  length  or 


greater  away  from  the  traction,  the  stress  for  four-point  integration 
has  an  error  of  less  than  0.6%  (Figure  6b).  At  this  distance  the  stress 
is  adequately  computed  by  four-point  quadrature.  The  region  closer  to 
the  element  is  herein  referred  to  as  the  "error  region."  The  size  of 
this  region  is  dependent  on  the  acceptable  error  level.  The  oscillation 
observed  in  the  response  has  been  referred  to  as  the  ripple  effect  in 
previous  reports  (Ref  15  and  8). 

There  are  different  approaches  to  improving  the  accuracy  in  the 
error  region.  Some  approaches  reduce  the  size  of  the  error  region, 
while  others  avoid  using  the  boundary  integration  to  calculate  responses 
in  the  error  region.  Increasing  the  order  of  numerical  integration 
reduces  the  weight  factors  A^  and  the  distance  between  individual  Gauss 
points,  thus  reducing  the  size  of  the  error  region.  This  has  been 
illustrated  in  Figures  7  and  8,  which  give  the  stress  and  error  fields 
using  8-  and  16-point  numerical  integration,  respectively.  The  error 
region  is  sigificantly  reduced.  Other  approaches  for  reducing  the  error 
region  include  element  subdivision  and  analytical  integration  of  extracted 
portions  of  the  integrand.  The  element  subdivision  approach  involves 
subdividing  the  element  and  performing  the  integrations  over  smaller 
segments  of  the  element.  Alternatively,  shape  functions  can  be  used  to 
interpolate  through  the  error  region  (bridge  the  region),  since  boundary 
values  for  points  exactly  on  the  boundary  are  obtained  accurately. 

When  accuracy  requirements  exceed  the  order  of  integration  that  is 
available  in  the  program  (i.e.,  when  the  response  point  is  within  the 
error  region  of  a  given  element),  one  of  the  previously  mentioned  remedies 
may  be  used.  Interpolation  through  the  error  region  can  give  good 
results,  but  it  can  also  constrain  the  choice  of  the  element  size  if  the 
analyst  is  also  anticipating  a  need  for  calculating  near-boundary  response 
points.  From  the  analyst's  perspective  it  would  be  more  appropriate  to 
base  the  selection  of  element  size  only  on  considerations  of  geometry 
and  anticipated  response  gradients  along  the  boundary.  In  this  regard, 
an  element  subdivision  approach  is  attractive  since  it  obviates  the  need 
for  apriori  consideration  of  near-boundary  response  point  calculations 
when  designing  the  boundary  discretization. 


The  variation  in  size  of  the  error  region  with  order  of  numerical 
integration  suggests  that  the  order  of  integration  be  determined  as  a 
function  of  normalized  distance  to  the  element  r;  where  r  is  the  ratio 
of  the  distance  between  the  field  point  and  the  element,  r,  to  the 
element  length,  M.  Other  factors  include  the  element  distortion, 
traction  distribution,  and  the  position  of  the  response  point  relative 
to  the  element. 

In  this  study  the  order  of  integration  is  adaptively  selected  for 
each  combination  of  element  and  field  point.  It  is  given  by  a  procedure 
that  includes  r  and  a  user-specified  error  parameter,  based  on  work 
originally  done  by  Lachat  and  Watson  (Ref  5).  The  actual  formula  used 
to  select  the  order  of  integration  (M)  is  a  modification  of  one  given  by 
Banerjee  and  Butterfield  (Ref  10). 

The  variation  of  M  given  by  Banerjee  and  Butterfield's  formula,  for 
three  different  (user-defined)  error  thresholds,  is  shown  in  Figure  9. 

For  each  error  threshold  the  order  of  integration  increases  as  r  decreases 
until  M  unexpectedly  decreases  very  near  the  element.  The  decrease  in 
order  near  the  element  is  a  mistake  that  can  cause  less  reliability  in 
the  near-boundary  region.  Some  implementations  will  "key"  on  M  to 
determine  if  element  subdivision  is  necessary;  this  formula  would  then 
result  in  no  element  subdivision  where  it  is  needed  most.  Other  strategies 
for  refining  the  results  in  the  near-boundary  region  would  not  be  affected 
by  the  behavior  of  Banerjee  and  Butterfield's  formula. 

A  revised  version  of  the  formula  is  given  below  in  a  Pascal-like 
outline. 

PROCEDURE  ORDERcalc 
CONSTANTS 

ZEROtolerauce:=  -2.0E-3; 

GAUSS INTmax : =  16; 

BEGIN  (*  the  ORDERcalc  procedure  *) 

LNERR0Rover8:=  LN(ERR0Rparameter/8) 

LNLover4R:=  LN(LENGTH*0 . 25/RADIUSmin) 

IF  LNLove r4R>ZER0T0L  THEN  LNLover4R:=  ZEROTOL 
ORDER :=  R0UND(0.5*ABS(LNERR0Rover8/LNLover4R-l .0))  +  1 
IF  ORDER>GAUSSINTmax  THEN  ORDER :=  GAUSSINTmax 
END  (*  the  ORDERcalc  procedure  *) 


The  variation  of  M  given  by  this  formula,  for  three  different  error 
parameters,  is  shown  in  Figure  10.  The  revision  gives  the  formula  an 
asymptotic  behavior  as  the  element  is  approached.  In  the  research 
software  developed  for  this  study  (Single  Element  Integration  Tester 
(SEIT)  and  BEM  Quadratic  Element  TEST  (BQTEST)),  the  order  of  integration 
could  vary  from  1  to  16.  The  formula  merits  further  investigation. 

Though  it  is  a  simple  formula  it  can  drastically  affect  the  efficiency 
and  accuracy  of  the  method.  The  SEIT  program  could  readily  be  used  to 
fit  a  curve  to  the  numerical  error  data. 

Practical  limitations  to  increasing  the  order  of  numerical  inte¬ 
gration  can  arise,  especially  in  the  context  of  the  direct  boundary 

element  method,  where  the  singularity  for  stress  integrations  is  more 

2 

severe  than  in  the  indirect  method  (1/r  versus  1/r  ). 

The  element  subdivision  algorithm  developed  here  recursively  sub¬ 
divides  the  element.  It  is  implemented  in  an  adaptive  manner  and  results 
in  a  concentration  of  subelements  near  the  response  point  of  interest. 

The  use  of  recursion  allows  a  general  algorithm  that  is  otherwise  difficult 
to  implement  for  curved  elements  when  based  on  an  alternative  iterative 
approach.  The  adaptive  concentration  of  subelements  in  the  proximity  of 
the  response  point  helps  minimize  the  number  of  integration  points.  The 
following  section  gives  an  overview  of  the  algorithm  followed  by  a  more 
detailed  pseudo-code  (an  outline  of  an  algorithm  written  in  a  combination 
of  English  and  computer  programming  language). 

RECURSIVE  QUADRATIC  ELEMENT 

Consider  the  isoparametric  quadratic  element  shown  in  Figure  11. 

Based  on  experiments  illustrated  in  the  previous  section  a  constant 

ratio  r  ,  can  be  established  that  bounds  the  element  error  region  as 
sub 

shown  by  the  dashed  line.  In  our  studies,  for  example,  for  order  16 
integration,  responses  have  less  than  0.1%  error  for  values  of  r  >  0.2. 

This  value  is  obtained  for  an  element  with  both  a  constant  traction 
distribution  and  a  constant  Jacobian.  The  criterion  that  the  radius 
must  satisfy  to  avoid  subdivision  of  the  element  is  given  by 


r  >  M  r  ,  (19) 

sub 

Element  subdivision  must  be  used  to  obtain  consistent  accuracy  for 

response  points  inside  of  this  radius. 

The  radius  r  can  be  estimated  by  calculating  the  distance  to  the 

nearest  node  point  when  the  interior  response  point  position,  x?,  is 

relatively  far  from  the  element  (e.g.,  r  >  M) .  When  close  to  the 

element  the  radius  r  may  be  considered  as  a  continuous  function  of  the 

local  coordinate  r)-  The  minimum  radius  to  the  element  can  then  be 

calculated  by  minimizing  the  square  of  the  radius  with  respect  to 

The  closest  point  on  the  element  is  indicated  by  local  coordinate 

Hq  in  Figure  11a  and  is  the  focal  point  for  the  initial  element  sub- 

c  r 

division.  A  "critical"  element  AT  is  centered  at  x^(Hq)  as  shown  in 

c  r 

Figure  lib  with  a  length  of  A£  given  by 

MCr(C)  =  (20) 

rsub 

Assuming  a  constant  Jacobian  over  the  subelement,  the  length  in  the 
local  coordinate  is  given  as 


Afl 


cr 


MCr(£) 

J(n0) 


(21) 


The  inaccuracy  introduced  by  the  assumption  of  a  constant  Jacobian  can 

be  overcome  if  the  user-specified  ratio  r  ^  is  set  slightly  greater 

than  the  value  computed  based  on  this  assumption  (e.g.,  0.3  instead  of 

cr 

0.2  in  this  case).  This  can  result  in  more  subdivision  since  A Z  will 


*The  possibility  of  obtaining  two  minimum  radius  values  poses  no 
problem  to  the  recursive  procedure. 
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be  smaller,  but  the  inclusion  of  variable-order  integration  will  reduce 
the  order  of  numerical  integration  on  the  individual  subelements  if  the 
specified  ratio  r  ^  is  overly  conservative. 

The  element  integrations  can  now  be  performed  with  consistent 

cr  cr 

accuracy  over  AT  with  respect  to  the  subelement  local  coordinate  q 

c  r 

and  subelement  shape  functions  N^(q  ).  The  contribution  of  subelement 
ATcr  in  element  q  to  the  traction  t?,  for  example,  at  x?  is  given  by 


M=3 
1  P 
a=l 


ka 


M=3 
I  N 


S  Nfl(ncr)  j  dncr 


Zj  ik  "p 


(22) 


where  N  _  is  the  value  the  oth  element  shape  function  at  node  B  of  the 

Ofp 

subelement.  The  quadratic  element  shape  functions  are  thus  being 

approximated  by  the  quadratic  subelement  shape  functions  N_.  In  a 

P 

similar  manner,  integrations  over  the  remainder  of  the  element  sub- 
3  b 

regions  AT  and  AT  shown  in  Figure  lib  can  be  incorporated  into  the 
original  element  integrations.  Expression  22  is  the  key  to  subdivision 
algorithms,  whether  they  are  implemented  with  iterative  or  recursive 
algorithms . 

The  minimum  radius  criterion,  Equation  19,  is  also  applied  to 
3  b  cr 

subregions  AT  and  AT  .  Subelement  AT  explicitly  satisfies  the  crite¬ 
rion  by  Equation  20.  The  remaining  element  subregions,  in  general,  do 
not  satisfy  the  subdivision  criterion  and  must  be  further  subdivided.  A 
recursive  application  of  the  subdivision  procedure  described  above  can 
then  be  applied  to  the  remaining  element  subregions.  Figure  11c  illus- 

trates  this  process  by  subdividing  AT  .  Primes  denote  one  level  of 

£ 

recursion.  The  closest  point  on  Ar  to  the  response  point  is  at  r]'  =  1. 

Since  this  is  at  the  end  of  the  subregion,  the  critical  subelement  is 

not  centered  on  q' ,  and  further  it  is  assumed  that  AT^  in  this  case  is 
”  b 

not  needed,  (i.e.,  no  further  subdivision  of  AT  is  required).  The 

cr  * 

critical  subelement  AT  can  now  be  integrated  and  incorporated  into 

£ 

the  integrals  for  the  subregion  AT  shape  functions  using  Expression  22. 
£  * 

Assume  now  that  AT  satisfies  Equation  19  (i.e.,  no  further  subdivision 

£ 

is  necessary  for  AT  ) .  It  can  then  be  integrated  and  incorporated  into 

£ 

the  subelement  integrations.  With  the  integrations  complete  on  AT  , 
Expression  22  incorporates  contributions  into  the  integration  over  AT. 
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The  integrations  over  the  element  AT  using  the  original  element 
shape  functions  have  been  completed  except  for  the  element  subregion 

AT  .  As  indicated  in  Figure  lib,  the  length  of  AT*5  is  less  than  that  of 

c  r  |)  cr 

AT  ,  while  the  radius  r  to  AT  is  greater  than  that  to  AT  ;  Equation  19 

is  satisfied.  Thus,  in  this  case  AT  would  require  no  further  subdivision 

Subregion  AF^  can  be  integrated  and  incorporated  into  the  integrals  for 

the  element  shape  functions  using  Expression  22. 

In  software  programming  parlance,  recursion  is  the  ability  of  a 
function  or  subroutine  to  call  itself.  One  strength  of  a  recursive 
implementation,  as  contrasted  with  an  iterative  implementation,  is  that 
the  element  subdivision  algorithm  can  directly  deal  with  the  case  when 
the  point  q^  is  internal  to  AT  (not  at  an  extrema).  This  case  can 
occur  with  curved  elements.  An  iterative  implementation  of  the  subdivi¬ 
sion  procedure  would  require  significant  "bookkeeping"  to  monitor  which 
portions  of  the  element  AT  had  been  integrated.  The  recursive  approach 
also  yields  compact  code.  Recursion,  though  a  common  feature  in  "modern 
computer  languages"  such  as  Pascal,  Modula  II,  and  C,  is  not  a  standard 
feature  of  FORTRAN.  Pascal  was  the  programming  language  used  in  this 
study. 

An  outline  of  the  recursive  element  subdivision  algorithm,  written 
in  pseudo-code,  is  given  below.  In  this  description  in  pseudo-code, 

Pascal  procedure  (analogous  to  a  FORTRAN  subroutine)  calls  of  the  element 
are  either  (1)  replaced  with  a  description  of  the  procedure's  actions  or 
(2)  explicitly  shown  with  a  listing  of  the  procedure  after  the  main 
element  procedure  QUAD.  (A  procedure  call  is  not  prefaced  by  the  CALL 
statement  as  in  FORTRAN.)  Formal  parameters  of  the  procedures  are  not 
shown  in  Pascal  syntax  and  are  only  included  if  they  clarify  the  recursive 
nature  of  the  algorithm.  Comments  that  explain  rather  than  replace  code 
are  enclosed  within  the  delimiters  (*  and  *); 

PROCEDURE  QUAD 

(INTEGRATIONS,  (*  the  element  integration  results  *) 

A£,  (*  element  length  estimate  *) 

x)  (*  array  of  element  nodal  coordinates  *) 
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PROCEDURE  QUADnonsing 

(INTEGRATIONS,  (*  the  element  integration  results  *) 
M,  (*  element  length  estimate  *) 

x,  (*  array  of  nodal  coordinates  *) 

r)  (*  minimum  radius  to  the  field  point  *) 


BEGIN  (*  the  QUADnonsing  procedure  *) 

Calculate  the  order  of  numerical  integration  (M)  using  the  ORDERcalc 
procedure  previously  defined 
Intialize  the  INTEGRATIONS  array  to  zero 
FOR  i=  1  TO  M  (*  each  integration  point  *) 

BEGIN  (*  the  integration  point  loop  *) 

For  each  of  the  element  integrations  (e.g.,  as  in  Equation  17) 
add  the  product  of  the  integrand  at  the  Gauss  point  with  the 
Gauss  weight  to  the  appropriate  INTEGRATIONS  array  value 
END  (*  the  integration  point  loop  *) 

END  (*  the  QUADnonsing  procedure  *) 


BEGIN  (*  the  QUAD  procedure  *) 

Calculate  the  minimum  radius  (r)  and  the  close  point  ( rj0 ) 

IF  r>A£*rsuk  (*  subdivision  criteria  of  Equation  19  *) 

THEN  (*  the  element  does  NOT  need  to  be  subdivided  *) 

QUADnonsing(INTEGRATIONS,A£,x,r)  (*  integrate  the  entire  element  *) 

ELSE  (*  the  element  must  be  subdivided  *) 

BEGIN 

Initialize  the  INTEGRATIONS  values  to  zero 

Determine  preliminary  values  for  subelement  integrations 
cr 

calculate  by  Equation  20 

c  r 

calculate  x  an  array  of  the  critical  subelement  nodal  elements 
QUADnonsing(INTEGRATIONSCr ,MCr,xCr,r)  (*  integrate  the  critical 

subelement  *) 
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Update  the  element  integrations  by  Expression  22 
c  r  cr 

IF  x  ( rj j  )*x(r)j)  (*  subelement  node  1  does  not  coincide  with 

element  node  1  *) 

THEN  (*  a  AT  subelement  is  necessary  *) 

BEGIN 

Interpolate  the  coordinates  of  x  (2) 

3 

Estimate  the  length  A£  using  two-point  integration 

QUAD  ( INTEGRATIONS3,  AJ2a,x3)  (*  a  recursive  call  *) 

Update  the  element  integrations  as  Expression  22 

END 
c  r*  c  it 

IF  x  (n3  )^x(n3)  (*  subelement  node  3  does  not  conincide  with 

element  node  3*  ) 

THEN  (*  a  Afb  subelement  is  necessary  *) 

BEGIN 

Interpolate  the  coordinates  of  xb(2) 

Estimate  the  length  A£b  using  two-point  integration 
QUAD(INTEGRATI0NSb,A2b,xb)  (*  a  recursive  call  *) 

Update  the  element  integrations  as  Expression  22 
END 

END  (*  the  element  subdivision  *) 

END  (*  the  QUAD  procedure  *) 


The  recursive  procedure  outlined  above  applies  a  variable-order 
integration  formula  at  every  level  of  recursion.  The  effect  of  the 
algorithm  is  to  concentrate  integration  points  on  the  boundary  near  the 
interior  response  point  x?.  In  this  algorithm  integration  points  further 
from  x?  will  increase  in  weight  and  spacing.  This  is  consistent  with 
the  idea  of  variable-order  integration  applied  around  the  problem  boundary 
(Ref  5),  but  here  it  is  applied  at  the  element  level. 

The  maximum  order  of  integration  available  in  the  program  is  used 
in  the  determination  of  the  error  region.  If  the  error  region  were  to 
have  been  based  on  single-point  integration,  the  recursive  algorithm 
would  use  a  "minimum''  number  of  integration  points  adaptively  positioned 
to  account  for  the  element  distortion.  Though  the  number  of  integration 
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points  would  be  minimized,  the  number  of  recursive  subdivisions  would  be 
maximized.  Efficiency  studies  comparing  these  two  factors  are  a  future 
goal. 

Each  level  of  recursion  has  a  corresponding  amount  of  computational 

overhead.  Subelement  lengths  and  nodal  positions,  minimum  distance  r, 

integration  order  M,  and  incorporation  of  the  subelement  integration 

results  into  the  original  integrals  (Expression  22)  are  all  calculations 

associated  with  a  recursive  subdivision.  Each  level  of  recursion  also 

requires  memory  to  accommodate  the  local  variables  and  pointers  to 

formal  parameters  passed  to  the  recursive  procedure  (subroutine).  For 

small  values  of  r  the  number  of  recursion  levels  required  is  small, 

even  for  response  points  very  close  to  the  element.  For  a  constant 

Jacobian  element  with  r  ,  =  0.25,  the  maximum  number  of  levels  of 

sub 

recursion  for  r  =  0.01  and  r  =  0.001  is  2  and  3,  respectively. 

With  the  concepts  of  the  recursive  procedure  explained,  the  actual 
accuracy  of  the  recursive  element  is  assessed  in  the  next  section. 

NUMERICAL  RESULTS 

Consistent  with  the  previous  study  of  numerical  integration  error, 
the  accuracy  of  the  recursive  element  is  first  demonstrated  for  a  single 
element  in  an  infinite  plane.  Figure  12a  presents  the  stress  (a^)  in 
the  problem  defined  by  Figure  5.  Compare  these  results  with  those  given 
by  Figures  6  through  8.  The  recursive  element  removes  the  singularity 
behavior  previously  noted  at  the  Gauss  points.  Figure  12b  gives  the 
percent  error  in  the  integration.  Notice  that  with  r  ^  =  0.15,  there 
is  a  small  amount  of  ripple  present  for  r  >  0.15,  the  region  in  which 
subdivision  is  precluded.  For  r  <  0.15,  the  element  adaptively  sub¬ 
divides  the  element  and  the  error  is  virtually  eliminated.  If  r  ^  is 
increased  to  0.30,  the  error  is  further  reduced  but  at  the  cost  of 
increased  computational  effort.  If  the  radius  (r)  is  simply  computed 
using  only  the  minimum  nodal  radius,  the  results,  shown  in  Figure  13, 
are  poor  in  the  error  region  except  near  the  nodes  where  the  calculation 
method  is  valid.  Between  nodes,  r  is  over  estimated  and  thus  the  sub¬ 
division  criterion  of  Equation  19  is  mistakenly  satisfied. 


The  recursive  element  is  next  demonstrated  on  two  finite  domain 
problems.  The  classic  problem  of  a  hole  in  a  finite  plate  is  demon¬ 
strated  first.  This  test  problem  illustrates  the  accuracy  that  can  be 
attained  with  the  recursive  element  near  the  boundary  even  with  high 
stress  gradients  present.  The  second  test  problem  is  a  square  plate 
subjected  to  uniform  tension.  This  theoretically  simple  problem 
illustrates  the  difficulty  the  indirect  BEM  formulation  has  near 
geometric  discontinuities.  Also  a  parameter  study  was  conducted  to 
illustrate  how  the  error  field  is  affected  by  a  collocation  point 
location  near  the  corner. 

Figure  14  defines  the  stress  concentration  problem  and  indicates 
the  region  in  which  the  internal  response  is  calculated.  Sixteen 
quadratic  elements  were  used  to  model  the  boundary  of  the  plate:  eight 
for  the  square  and  eight  for  the  circle.  Discontinuous  elements  were 
used  in  the  corners.  The  collocation  points  for  the  "corner  nodes"  only 
were  positioned  at  |  H  |  =  0.80.  All  other  nodes  that  were  common  to  two 
elements  were  continuous  (i.e.,  the  collocation  points  were  placed  at 
the  geometric  nodes). 

In  developing  the  system  of  equations,  integrations  over  the  singu¬ 
larities  were  avoided.  The  collocation  points  were  positioned  at  a 
distance  r  =  0.001A£  from  the  element  inside  the  domain  ft.  The  mathe¬ 
matical  limit  was  thereby  replaced  by  a  physical  limit  since  the 
recursive  element  is  designed  to  handle  near-boundary  responses.  This 
is  not  necessarily  the  most  computationally  efficient  means  of  calculating 
the  integrations,  but  it  is  simple  to  implement  and  prevents  further 
complication  of  the  element  integration  procedure. 

The  calculation  of  internal  responses  is  performed  both  with  the 
simple  four-point  quadrature  method  and  with  the  recursive  technique. 

The  response  points  are  spaced  at  0.06  in  a  square  grid.  Thus,  the 
distance  of  the  response  points  from  the  circular  hole  varies  around  the 
circumference.  Figure  15  presents  the  stress  (°jj)  near  the  hole  using 
four-point  quadrature.  The  ripple  effect  near  the  hole  is  quite  prevalent 
The  same  responses  calculated  with  the  recursive  technique  are  given  in 
Figure  16.  The  oscillation  of  the  response  near  the  hole  is  removed; 
thus,  values  in  the  high  gradient  regions  of  interest  are  accurately 


predicted.  At  a  distance  of  0.001  from  the  top  of  the  hole  the  stress 
concentration  factor  was  3.06.  This  compares  well  with  theoretical 
solutions  for  an  infinite  plate,  3.00,  and  the  finite  plate,  3.04 
(Ref  16).  Error  with  respect  to  the  finite  plate  solution  (at  the  exact 
top  of  the  hole)  is  0.66%. 

Figure  17  defines  the  uniaxial  tension  problem  and  indicates  the 
region  in  which  the  internal  response  is  calculated.  The  problem  was 
modeled  using  both  four  and  eight  elements  of  uniform  length.  The 
collocation  positions  near  the  corner  were  varied  to  study  the  effect  on 
the  corner  response.  All  collocation  points  were  positioned  within  the 
domain  to  develop  the  system  of  equations  (as  explained  in  the  previous 
test  problem) . 

The  results  for  the  four-  and  eight-element  models  are  shown  in 
Figures  18  and  19,  respectively.  The  nonrecursive  integration  results 
use  the  variable-order  integration  formula  but  are  limited  to  a  maximum 
integration  order  of  16.  As  with  the  previous  problem  the  recursive 
element  effectively  eliminates  the  near-boundary  error  that  is  charac¬ 
teristic  of  the  usual  implementation  of  Gauss  quadrature.  The  computed 
artificial  boundary  tractions  for  both  models  are  sketched  in  Figure  20. 
The  improvement  in  modeling  these  artificial  tractions  with  an  increase 
in  the  number  of  elements  was  apparent  in  the  computed  stress  response 
of  Figures  18b  and  19b.  In  the  previous  test  problem,  responses  were 
examined  near  a  geometrically  continuous  portion  of  the  boundary  (the 
hole)  so  no  severe  gradients  occur  in  the  corresponding  artificial 
tractions.  However,  in  this  test  problem  the  responses  approach  the 
corner  that  has  geometry  and  traction  discontinuities  and  gradients  are 
present  in  the  corresponding  artificial  tractions.  The  error  in  the 
corner  reflects  an  inherent  weakness  of  the  indirect  boundary  element 
formulation.  The  artificial  boundary  tractions  become  unbounded  near 
the  corner  and  are  difficult  to  represent  with  polynomial  element  shape 
functions . 

The  positions  of  the  corner  collocation  points  and  the  element 
middle  nodes  (eight-element  model)  can  be  varied  to  improve  the  results 
without  significantly  affecting  the  cost  of  the  analysis.  The  effect  of 
shifting  the  middle  nodes  toward  the  corner  would  probably  improve  the 


results  since  the  middle  node  collocation  point  would  be  in  a  position 
to  better  model  the  gradient.  The  middle  node  position  was  not  considered 
in  this  study.  By  similar  reasoning,  the  position  of  the  discontinuous 
corner  node's  collocation  point  is  expected  to  affect  the  results. 

Figure  21  gives  the  error  field  in  the  domain  for  various  positions  of 
the  corner  collocation  points.  The  error  field  for  this  problem  is 
minimized  when  r)  0.80.  The  error  approaches  zero  near  the  colloca¬ 
tion  points  where  the  boundary  conditions  are  satisfied  exactly.  In 
addition  to  the  difficulty  of  modeling  the  artificial  tractions,  col¬ 
location  points  positioned  very  close  to  the  corner  can  result  in 
equations  that  approach  linear  dependence  and  are  numerically  ill- 
conditioned. 

CONCLUSIONS 

A  new  recursive  element  subdivision  algorithm  is  presented  for  the 
boundary  element  method.  The  algorithm  is  introduced  within  the  context 
of  the  indirect  boundary  element  method  using  an  isoparametric  quadratic 
element;  it  is,  however,  generally  applicable  to  the  direct  boundary 
element  method  and  other  element  formulations. 

The  motivation  for  this  study  and  the  resulting  recursive  algorithm 
development  is  to  improve  the  reliability  of  the  boundary  e1ement  method 
by  improving  accuracy  in  the  near-boundary  region  where  error  is  normally 
excessive.  As  a  result  of  the  algorithm,  the  analyst  can  choose  element 
size  on  the  basis  of  geometry  and  anticipated  response  gradients  alone 
without  concern  for  avoiding  the  error  in  the  internal  response  in  the 
near-boundary  region.  This  is  similar  to  the  idea  behind  the  isopara¬ 
metric  element  formulations  in  the  finite  element  method,  where  fewer 
elements  are  necessary  to  capture  geometric  description  while  at  the 
same  time  maintaining  sufficient  accuracy  to  model  higher  gradients. 

In  the  indirect  boundary  element  formulation  the  gradients  of  the 
artificial  boundary  tractions  must  be  considered  when  subdividing  the 
boundary.  This  is  a  practical  limitation  of  the  indirect  method  since 
these  gradients  are  more  difficult  to  identify  on  a  physical  basis.  As 


illustrated  in  this  study,  the  artificial  tractions  can  have  very  high 
gradients  near  geometric  and  loading  discontinuities.  Concentrating 
elements  near  the  discontinuities  is  one  way  to  capture  these  gradients. 

The  computational  liability  attending  the  new  algorithm  is  the 
recursive  subdivision  of  a  few  elements,  but  this  is  at  least  partially 
offset  by  a  reduction  of  elements  otherwise  required  in  anticipation  of 
near-boundary  response  point  calculations. 

The  recursive  algorithm  is,  on  the  basis  of  preliminary  implementa¬ 
tion,  believed  to  be  relatively  simple  to  implement  compared  to  an 
iterative  approach  that  would  also  allow  the  same  generality  with  curved 
elements.  It  also  lends  itself  to  surface  elements  for  three-dimensional 
applications . 

The  adaptive  quality  of  the  recursive  subdivision  procedure  poten¬ 
tially  enhances  the  use  of  the  method  in  computer-aided  design  (CAD) 
environments  where  the  user  may  not  have  a  complete  understanding  of  the 
numerical  behavior  of  the  boundary  element  method  in  the  near-boundary 
region.  It  would  appear  that  it  is  also  applicable  to  stress  analysis 
applications  in  plasticity  and  fracture  mechanics  where  response  accuracy 
near  the  boundary  is  vitally  important. 

RECOMMENDATIONS 

One  of  the  objectives  in  our  research  has  been  a  combined  finite 
element  and  boundary  element  program  that  could  reduce  the  high  cost  now 
associated  with  nonlinear  finite  element  programs.  Further  work  is 
needed  for  an  effective  coupling  with  nonlinear  finite  element  programs. 

A  boundary  element  formulation  that  produces  a  symmetric  set  of  equations 
is  essential.  The  concepts  for  an  approach  that  produces  a  symmetric 
stiffness  matrix  were  outlined  in  this  year’s  effort.  This  approach 
would  allow  a  much  easier  coupling  of  the  BEM  with  existing  FEM  programs. 

The  recursive  algorithm  in  this  year's  effort  adaptively  integrated 
the  elements.  The  adaptive  refinement  of  a  combined  FEM/BEM  "mesh" 
could  drastically  increase  the  efficiency  of  this  class  of  problems. 

The  application  of  a  recursive  algorithm  to  adaptively  refine  a  FEM/BEM 
mesh  should  be  investigated. 


The  authors  recommend  that  two  areas  receive  further  investigation: 
(1)  a  stiffness  coupling  of  the  BEM  with  the  FEM,  and  (2)  adaptive 
methods  for  combined  boundary  and  finite  elements. 
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(a)  stress  field. 


(b)  CTji  error  field. 

Figure  7.  Eight-point  integration  results  for  a  single  element 
in  an  infinite  plane  (defined  in  Figure  5). 
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error  region 


(a)  Critical  point  for  element  subdivision. 


Figure  13.  Recursive  element  subdivision. 
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Figure  17.  Square  plate  in  uniaxial  tension  problem. 
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(a)  CTji  error  field  without  recursive  integration. 


4 


42 


INSTRUCTIONS 


The  Naval  Gvil  Engineering  Laboratory  has  revised  its  primary  distribution  lists.  The  bottom  of 
the  mailing  label  has  several  numbers  listed.  These  numbers  correspond  to  numbers  assigned  to  the  list  of 
Subject  Categories.  Numbers  on  the  label  corresponding  to  those  on  the  list  indicate  the  subject  category  and 
type  of  documents  you  are  presently  receiving.  If  you  are  satisfied,  throw  this  card  away  (or  file  it  for  later 
reference). 

If  you  want  to  change  what  you  are  presently  receiving. 

•  Delete  —  mark  off  number  on  bottom  of  label. 

•  Add  -  circle  number  on  list. 

•  Remove  my  name  from  all  your  lists  -  check  box  on  list. 

•  Change  my  address  -  line  out  incorrect  line  and  write  in  correction  (ATTACH  MAILING  LABEL). 

•  Number  of  copies  should  be  entered  after  the  title  of  the  subject  categories  vou  select. 

Fold  on  line  below  and  drop  in  the  mail. 

Note:  Numbers  on  label  but  not  listed  on  questionnaire  are  for  NCEL  use  only,  please  ignore  them. 


Fold  on  Ime  and  staple. 


DEPARTMENT  OF  THE  NAVY 


NAVAL  CIVIL  ENGINEERING  LABORATORY 
PORT  HUENEME.  CALIFORNIA  93043 


POSTAOE  AND  FEU  PAID 
DEPARTMENT  OF  THE  NAVY 
DOD-Sta 


OFFICIAL  BUSINESS 
PENALTY  FOR  PRIVATE  USE.  SSOO 
I  IND.NCEL.2700/ 4  (REV.  12.73) 
0SS0-LL-L7O.0044 


Commanding  Officer 
Code  LI  4 

Naval  Civil  Engineering  Laboratory 
Port  Hueneme,  California  93043 


DISTRIBUTION  QUESTIONNAIRE 

The  Naval  Civil  Engineering  Laboratory  is  revising  its  primary  distribution  lists. 


SUBJECT  CATEGORIES 

1  SHORE  FACILITIES 

2  Construction  methods  and  materials  (including  corrosion 

control,  coatings) 

3  Waterfront  structures  (maintenance/deterioration  control) 

4  Utilities  (including  power  conditioning) 

5  Explosives  safety 

6  Construction  equipment  and  machinery 

7  Fire  prevention  and  control 

8  Antenna  technology 

9  Structural  analysis  and  design  (including  numerical  and 

computer  techniques) 

10  Protective  construction  (including  hardened  shelters, 

shock  and  vibration  studies) 

1 1  Soti/rock  mechanics 

13  BEQ 

14  Airfields  and  pavements 

15  ADVANCED  BASE  AND  AMPHIBIOUS  FACILITIES 

16  Bate  facilities  (including  shelters,  power  generation,  water  supplies) 

17  Expedient  roads/a irfieids/bndges 

18  Amphibious  operations  (including  breakwaters,  wave  forces) 

19  Over  the  Beach  operations  (including  containerization, 

materiel  transfer,  lighterage  and  cranes) 

20  POL  storage,  transfer  and  distribution 

24  POLAR  ENGINEERING 

24  Same  as  Advanced  8ase  and  Amphibious  Facilit  es, 
except  limited  to  cold  region  environments 


TYPES  OF  DOCUMENTS 

85  Tech  data  Sheets  86  Technical  Reports  and  Technical  Notes 

83  Table  of  Contents  &  Index  to  TDS 


28  ENERGY/POWER  GENERATION 

29  Thermal  conservation  (thermal  engineering  of  buildings,  HVAC 

systems,  energy  loss  measurement,  power  generation) 

30  Controls  and  electrical  conservation  (electrical  systems, 

energy  monitoring  and  control  systems) 

31  Fuel  flexibility  (liquid  fuels,  coal  utilization,  energy 

from  solid  waste) 

32  Alternate  energy  source  (geothermal  power,  photovoltaic 

power  systems,  solar  systems,  wind  systems,  energy  storage 
systems) 

33  Site  data  and  systems  integration  (energy  resource  date,  energy 

consumption  data,  integrating  energy  systems) 

34  ENVIRONMENTAL  PROTECTION 

35  Solid  waste  management 

36  Hazardous/toxic  materials  management 

37  Wastewater  management  and  sanitary  engineering 

38  Oil  pollution  removal  and  recovery 

39  Air  pollution 

40  Noise  abatement 

44  OCEAN  ENGINEERING 

45  Seafloor  soils  and  foundations 

46  Seafloor  construction  systems  and  operations  (including 

diver  and  manipulator  tools) 

47  Undersea  structures  and  materials 

48  Anchors  and  moorings 

49  Undersea  power  systems,  electromechanical  cables, 

and  connectors 

50  Pressure  vessel  facilities 

51  Physical  environment  (including  site  surveying) 

52  Ocean-based  concrete  structures 

53  Hyperbaric  chambers 

54  Undersea  cable  dynamics 


82  NCEL  Guide  8c  Updates 
91  Physical  Security 


□  None- 

remove  my  name 


ov 


